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The peeling theorem of general relativity predicts that the Weyl curvature scalars (n = 
0, . . . ,4), when constructed from a suitable null tetrad in an asymptotically flat spacetime, fall off 
asymptotically as r""'' along outgoing radial null geodesies. This leads to the interpretation of *I'4 
as outgoing gravitational radiation at large distances from the source. We have performed numerical 
simulations in full general relativity of a binary black hole inspiral and merger, and have computed 
the Weyl scalars in the standard tetrad used in numerical relativity. In contrast with previous 
results nn , we observe that all the Weyl scalars fall off according to the predictions of the theorem. 

PACS numbers: 04.25.dg, 04.20.Ha, 0430.-W 



I. INTRODUCTION 

Spacetimes containing multiple black holes (BHs) are 
perhaps the most studied class of non-perturbative, low- 
symmetry vacuum solutions of Einstein's equations and 
have already opened the path to tests of the generality 
of a number of fundamental theorems and conjectures in 
general relativity. 

Studies of multiple-BH systems have addressed the be- 
haviour of trapped surfaces and tubes [2H8] and event 
horizons |9Hl5j during the dynamical many-body phase. 
The asymptotic properties (such as the radial falloff of 
the Weyl scalars) of binary-black-hole (BBH) spacetimes 
have been investigated p] , and the algebraic character of 
BBH spacetimes has been studied 16, 17J. A program 
to include conformal infinity in the numerical domain 
has been carried out [HI [Hj , and properties of higher- 
dimensional BBH systems have been described . 20-ri22] . 
Several interesting directions remain to be explored in 
the context of multiple-BH evolutions [^51128[ . 

It seems fair to state that three-dimensional numeri- 
cal evolutions of Einstein's equations are now leading to 
compelling and useful results for gravitational fields in 
non-trivial configurations. Whilst this class of numerical 
relativity (NR) studies may appear somewhat removed 
from the field of gravitational- wave astrophysics, the clar- 
ification of fundamental issues has a direct impact there. 
For instance, the interpretation of as the outgoing 
component of gravitational radiation relies entirely on 
the applicability of the peeling theorem. 

In this work, we follow up the study carried out in [T] , 
where the radial falloff of the Weyl scalars was mea- 
sured in a BBH spacetime and compared with the re- 
sults expected from the peeling theorem. We use the 
same evolution code and the same definitions of the Weyl 
scalars. However, we use an independent and open anal- 
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ysis framework provided by the Einstein Toolkit initia- 
tive [29 for computing the scalartj^ and we analyze their 
falloff properties using a more direct approach. Our re- 
sults indicate that all of the fall off at the expected 
rates, suggesting that the usual techniques used for wave 
extraction in NR simulations are sufficient for recover- 
ing the results of the peeling theorem. We present our 
method in Sec. HI] and the results we obtain in Sec. IIIII 
We conclude with a discussion in Sec. |IV[ Throughout 
this paper, we use a spacelike signature (— , +, +, +) and 
a system of units in which c = G = 1. 

II. METHODOLOGY 

A. The Weyl scalars and the peeling theorem 

The work of Sachs [31], Newman and Penrose [32] il- 
lustrates how the information encoded in the traceless 
part of the Riemann curvature tensor, the Weyl tensor 
Cabcd, can be completely expressed as a set of five com- 
plex numbers usually referred to as the Weyl scalars: 

*o - Cabcdrm'>tm'' (la) 

*i = Cabcd''n''e'm'' (lb) 

*2 = C.b.drm^m^n'^ (Ic) 

*3 - Cabcde''n''m'n'' (Id) 

^4 = Cabcdn''m''n'=m'' (le) 

where {n°' , , m°' , ffi"-) is a tetrad of two real and two 
complex null vectors satisfying nau"- = 0, (af"" = 0, 
Uai^ — —1, rriam"' = and mam"' = 1. 



^ The unexpected results pointed out in [T] for the peeling proper- 
ties of 'I'o and 'ti were affected by errors in the original analysis 
implementation which have subsequently been corrected. As an 
aside, we note that ensuring correctness is a subtle and involved 
part of developing complex scientific applications 1301 . We be- 
lieve that collaborative development and the use of open-source 
software is one of the simplest and most effective ways to ensure 
code correctness. 
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Handling the curvature in this notation has a number 
of benefits: from the geometric standpoint, there is im- 
mediate insight to be gained as the five are simply the 
spin-frame components of the Weyl spinor "^abcd- Thus 
if the spin-frame has some physical relevance (e.g. is ori- 
ented along principal null directions), then the can 
be identified as specific (radiative vs. Coulomb, trans- 
verse vs. longitudinal) components of the gravitational 
field [33]. From the practical side, the five complex quan- 
tities above simplify the analysis of the asymptotic prop- 
erties of the Ricmann tensor (see, e.g., [321[33]). In partic- 
ular, the work of Newman and Penrose shows how these 
scalars fall off along outgoing radial null geodesies in a 
neighbourhood of future null infinity, J^^, as: 



^ r. 



(2) 



The usual assumptions of the Penrose conformal con- 
struction are assumed, namely the existence of a space- 
time conformally related to the physical one by g^^^^^^^ = 
^^9ab^^' with rt — a.t ^ the falloff of the correspond- 
ing Weyl tensor at least as C^^^^ ^ il, also at , and 
the alignment of n° and 1°" with the ingoing and outgoing 
null directions, respectively (accomplished, for instance, 
by setting Ua — ri;a)- In the above, r = fi^^ can be used, 
at least in a neighbourhood of , as the affine param- 
eter on the null geodesic. In typical asymptotically-flat 
spacetimes, this can be identified with a radial coordi- 
nate, as long as the Weyl tensor satisfies the above falloff 
condition when expressed in terms of it. This result, 
known as the peeling theorem, is a convenient tool for 
evaluating integrals at J^~^ and allows one to identify ^'4 
with the outgoing gravitational-wave degrees of freedom 
since, for one, it is the only component with an associated 
fiux at J""*" that is not identically zero: 



4. 



(3) 



J+ 



Notice that, whilst this result is rather robust and valid 
in a large class of tetrad (or spinor) frames, it does rely 
on a suitable choice of r. A question that was recently 
raised in ^1] is to what extent the peeling theorem can be 
directly applied to the Weyl scalars usually calculated in 
NR simulations of BBH systems, where a dynamically- 
evolving gauge can obscure the character of the coordi- 
nates, from which the tetrad is usually constructed. This 
question is closely related to similar investigations carried 
out in the NR literature, such as what conditions should 
be imposed on the tetrad frame in order to retain the in- 
terpretation of \I'4 as outgoing gravitational waves [57], or 
how many principal null directions exist in post-merger 
black hole remnants [Ml [17] . 

In this work, we investigate the extent to which the 
falloff rates predicted by the peeling theorem apply to 
BBH spacetimes as usually computed in NR. There are 
two aspects to this question: whether a numerical space- 
time satisfies the theorem's physical requirements (the 
falloff condition for the Weyl tensor), and whether the 



dynamical coordinates are compatible with the assump- 
tions. 



B. Tetrad 

In numerical simulations, where the domain does not 
typically include J^~^ (the notable exception being simu- 
lations using characteristic extraction [Ullin]), the tetrad 
(n", m°, m") is usually chosen according to the follow- 
ing straightforward prescription [35 . First, one defines 
the spatial vectors 0", and 6*° via their spatial com- 
ponents: 



= [x,y,z], 
0' = V7e%-fc'/'^>-'- 



(4a) 
(4b) 
(4c) 



where 7 is the determinant of the spatial metric. One 
then uses Gram-Schmidt orthonormalisation to produce 
the orthonormal triad 



II' \\r-Pe,rr 



P.JW 



(5) 



where PxH"^ = {xi,y'')x'^ is the projection of a vector 
along a unit vector a;°, ||u|| — ^/vaV^, and indices are 
raised and lowered with the spatial metric. Finally, one 
complements this triad with the unit hypersurface nor- 
mal and constructs the null tetrad as 



1 



1 



V2' V2 
m" = 4=(e^ + *e;), m"^ = -^(eg - ze^). (6) 



V2' 



The normalisation in Eq. ^ follows a common conven- 
tion [3S] used in NR and differs from the conventions usu- 
ally used in other fields. This can lead to constant-factor 
differences between expressions computed using this con- 
vention and those appearing elsewhere in the literature. 

Unfortunately, due to the degeneracy in the azimuthal 
coordinate, 0, this tetrad is not defined on the z— axis. 
We overcome this problem by making the particular 
choice (/) = TT, so that 0* = [1, 0, 0] and the tetrad is once 
again well defined (note, however, that the tetrad — and 
hence the Weyl scalars — remain discontinuous across the 
z-axis) . 



C. Approximate null geodesies 

We wish to investigate the extent to which the Weyl 
scalars as computed using the tetrad in Sec. jll Bj obey the 
peeling property along some suitably defined curves in a 
BBH spacetime. For the peeling theorem to apply, the 
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curves should be outgoing radial null geodesies; this is 
because the Weyl scalars vary not only in advanced time 
(according to the peeling theorem), but also in retarded 
time. A failure to use exact outgoing null geodesies as 
the curves along which to measure the falloff results in 
a mixing of these two dependencies, leading to an error 
in the measured falloff. This error may be significant 
or not, depending on the extent of the deviation from 
exact geodesies. It would be possible to calculate the ex- 
act null geodesies in a numerical spacetime, but this is 
not typically done, and we have not done so here as it 
turns out to be unnecessary for confirming the applicabil- 
ity of the peeling theorem to BBH spacetimes. Instead, 
we compute a series of approximate null geodesies, de- 
fined as curves in the x-t coordinate plane along which 
local maxima of 1^*4 1 propagate. In the geometric optics 
approximation, these curves will be null geodesies. 

By restricting the approximate null geodesies to the x— 
t plane, we are assuming that the outgoing null geodesies 
which are asymptotically radial have a negligible angular 
coordinate dependence in the region in which we study 
them. This is a reasonable expectation given that in 
Kerr, with spin parameter equal to that of the final post- 
merger black hole {J/IvP ~ 0.68), the angular deviation 
of an asymptotically radial outgoing null geodesic be- 
tween r — BOM and r — lOOOM is only 10"'^ radians. 

One may expect that an even simpler approxima- 
tion would suffice and make the assumption that the 
null geodesies correspond approximately to those of the 
Schwarzschild spacetime in Schwarzschild coordinates. 
In this case the null geodesies are given by = t-f- const, 
where r* = r -I- 2Mln(r/2M — 1) is the tortoise coordi- 
nate. At large radius, this is not an unreasonable expec- 
tation. However, close to the black holes this approxima- 
tion becomes increasingly poor. Unfortunately, ^'o and 
'i'l fall off sufficiently fast that they drop below the level 
of numerical error before the region r « lOOM where the 
Schwarzschild approximation becomes reliable, so this 
approximation is unsuitable for measuring the falloff of 
these scalars. 



D. Falloff measurement 

Given the Weyl scalars computed using the coordinate 
tetrad ^ and an approximate null geodesic r = A(i), 
described in the previous section, we consider the falloff 
of the scalars with the coordinate radius r. In order to 
determine a representative value for the falloff rate, we 
choose an interval r e [rmin , ''max] and perform a linear 
least-squares fit of 

log|*„(i,A(t))| = -plogr-^ const (7) 

to determine p. 

The appropriate choice of the interval is influenced by 
two factors. Firstly, since the peeling theorem predicts 
only the asymptotic behaviour of the scalars, we desire a 
fitting interval at as large a radius as possible so that the 



mh/M mb/M D/M pr/M pt/M 
"as 0.47653463302 6 -0.005867766 0.138357448 



TABLE 1. Initial data parameters for the equal-mass, non- 
spinning configuration studied. M is the sum of the initial 
irreducible masses of the black holes, and mh is the irreducible 
mass of each BH. mb/M, D/M, pr/M and pt/M are the bare 
mass, separation, radial and tangential momenta used in the 
Bowen-York initial data prescription. 

sub-leading terms are negligible. Secondly, we find that 
^'o and ^'i are dominated by numerical error at large r 
and large t, and that this limits the maximum radius at 
which we can fit. 



E. Numerical implementation 

Our numerical BBH solutions were obtained using 
the Llcuna code [36] to solve Einstein's equations in 
the BSSN [37H39] formalism with the moving punc- 
ture method [101 SI] using eighth-order finite differenc- 
ing. The computational infrastructure is based on the 
Cactus framework HU |33] and the Carpet [ilHiS] adap- 
tive mesh-refinement driver, and implements a system of 
multiple grid patches with data exchanged via interpola- 
tion [IJB^. This multipatch technique allowed us to use a 
spherical outer grid with constant angular resolution to 
best match the resolution requirements of radially out- 
going waves, leading to an outer boundary at very large 
radius at only modest computational cost. This enabled 
us to measure the Weyl scalars accurately at a large ra- 
dius, where they are closer to their asymptotic form. 

We studied an equal-mass binary system with an ini- 
tial separation of 6M and both black holes initially non- 
spinning. The configuration studied in [Tj had a larger 
initial separation and hence a longer gravitational wave 
signal. However, the length of the wave signal should 
not be important for measuring the falloff, so we ran a 
shorter simulation for computational efficiency. We com- 
puted standard Bowen-York jlTl HH] initial data with pa- 
rameters given in Table |l] using the TwoPunctures [4^ 
code. 

In our simulations, calculation of the Weyl scalars 
was performed by WeylScal4, a component of the open 
Einstein Toolkit |5D]. We have adapted WeylScal4 
to utilise the Llama multipatch framework and have 
performed careful testing against analytic results for 
the /-invariant |35j in Kerr, as well as cross-checking 
WeylScal4 against the Psikadelia code. 

We performed simulations at two different grid spac- 
ings (the detailed grid structures used in our runs are 
listed in Table |ll]) to assess the effect of the numeri- 
cal discretisation on the solution. Owing to high fre- 
quency numerical refiections from the inter-patch and re- 
finement boundaries, we found it necessary to modify the 
grid structure when computing the falloff of and ^'i 
which, being the fastest decaying and lowest in ampli- 
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Kun 


ho/M Wang 


Hin/M 


Hout/M 


A'lov 


ri/M 


Ma 


0.69 28 


39.77 


2697.6 


6 


12,6,3,1.5,0.6 


Mb 


0.60 32 


39.6 


2697.6 


6 


12,6,3,1.5,0.6 


Ca 


0.69 




200.00 


6 


12,6,3,1.5,0.6 


Cb 


0.60 




200.00 


6 


12,6,3,1.5,0.6 



TABLE II. Numerical grid parameters of the BBH simulations 
studied, ho is the grid spacing on the coarsest Cartesian grid. 
In cases where an angular grid is used this is also equal to the 
radial grid spacing in the angular patches. A^'ang is the num- 
ber of cells in the angular directions in the angular patches. 
i?in and -Rout are the inner and outer radii of the angular 
patches. Ni^v is the number of refinement levels (including 
the coarsest) on the Cartesian grid, and ri indicates that a 
cubical refinement box of side 2ri is centred on the BH on 
level "1", with level being the coarsest. 



tude, are the most affected by these reflections. In this 
case, we eliminated the angular grid patches entirely and 
used a purely Cartesian domain with an outer bound- 
ary causally disconnected from the region in which we 
could resolve these Weyl scalars. We found that the er- 
rors coming from refinement boundaries visible in 
could be reduced significantly by switching the mesh 
refinement algorithm used from the standard Berger- 
Oliger [5T] scheme to the tapered grids scheme of [5^ . 
This has the effect of eliminating the errors caused by 
time interpolation at the mesh-refinement boundaries. 



III. RESULTS 

We now present the results of applying the fitting 
methods to the simulations described in Sec. Ill El 

We consider the falloff of the Weyl scalars along curves 
Xi in spacetime which correspond to the peaks i — 1 . . .7 
of |vl'4| shown in Fig. [I] This is a plot of |5'4| along the 
X-axis at fixed coordinate time t = 720M. Note that 
for this configuration, all of the Weyl scalars are either 
purely real or purely imaginary along the x-axis. 

Figure [2] plotted on a log- log scale, shows the falloff of 
the Weyl scalars along the curve Ai in the x-t plane cor- 
responding to the peak labelled 1 in Fig. [T] as a function 
of the radial coordinate r. The peeling theorem predicts 
that the Weyl scalars ^„ fall off asymptotically as r"~^ 
corresponding to straight lines in this figure. The solid 
curves represent the numerical data and the dashed lines 
represent the best fit straight line (the line for the ex- 
pected falloff is not shown, but in each case is visually 
very similar to the line obtained from the fit). We see 
that to a very good approximation, each scalar exhibits 
a power-law decay as predicted by the peeling theorem. 

For each scalar, we choose an appropriate fitting inter- 
val (see Sec. II D) and compute the falloff rate by per- 
forming a least-squares fit in this interval. For '^2, '^3 
and 4'4, the scalars are well resolved at all radii in our 
simulations (we stopped the simulations when the sig- 
nal reached r = lOOOM). We therefore choose a fitting 




FIG. 1. 1^4! on the i-axis. Note that since there is no mul- 
tipolar decomposition, this looks different to the waveform 
plots usually shown in papers. Specifically, the junk radia- 
tion is on the right and the wave propagates to the right. 
The numbered labels indicate the peaks corresponding to the 
approximate null geodesies Xi that we compute. 



interval r e [lOOM, lOOOM] in this case. For ^0 and 
^1, the scalars are dominated by numerical error beyond 
r « 60M, as we observed when comparing the results 
at different numerical resolutions. We therefore choose 
a fitting interval r G [30M,45M] for these. The fitting 
intervals are indicated on the plot as black vertical lines 
on each curve. 

We have performed this analysis for each of the 7 
curves and for each resolution. For ^'o and only Ai 
and A2 yield meaningful measurements of the falloff rate; 
there is too much finite differencing error in the solution 
on A3-A7 to compute a falloff rate. 

The shaded regions in Fig. [2] (visible only for ^ffQ and 
^1) represent an indication of the error due to finite dif- 
ferencing, computed using the two different resolutions 
and assuming eighth-order convergence of the error. This 
should be taken as indicative only, as the convergence 
rate obtained can vary between 3 and 8 at finite resolu- 
tion, corresponding to the different sources of numerical 
error in the simulation. 



Table III shows the falloff rates for each scalar along 
Ai as well as the rate quoted in [1] and the rate expected 
from the peeling theorem. The finite differencing error in 
the last digit is indicated in parentheses and again should 
be taken as only a coarse estimate. The rates obtained 
from A2-A7 (not shown) do not differ significantly from 
those obtained for Ai. Our measured rates are within 4% 
of the values expected from the peeling theorem. 

By studying the peeling properties of each of the Weyl 
scalars we may gain insight into where 4*4 may be used 
as a reasonable measure of the gravitational wave sig- 
nal. This is closely related to the identification of the 
regions referred to as near zone, transition zones and ra- 
diation zone in |32j (note that these are not the same 
zones referred to in [SB]). We illustrate this visually (for 
the curve for Ai) in Fig. [3j where we plot the relative 
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Expected 


Falloff rate p 
Ref. Oj 


Measured 





5 


2.00 


4.8(4) 


1 


4 


2.48 


3.91(7) 


2 


3 


2.99 


2.99307(6) 


3 


2 


1.99 


2.0135(6) 


4 


1 


0.99 


1.01333(7) 



TABLE III. FallofF rates for the Weyl scalars, including the 
rate expected from the peehng theorem, the rate obtained 
in [l], and the rate measured from our simulations for the 
approximate null geodesic Ai. The error in the last digit in- 
dicated in parentheses is a coarse estimate of the finite differ- 
encing error in the falloff rate. 
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FIG. 2. The falloff of the Weyl scalars along the approximate 
null geodesic Ai corresponding to the first inspiral peak of 
on the X-axis. The falloff rates obtained from fitting in 
the intervals indicated by vertical black lines are indicated in 
Table Hn] 



contribution from each of the Weyl scalars to the total 
curvature^ which we define as \^ = E„=ol*nl- Each 
shaded region Z„ indicates the region in which all the 
^'fc with k > n give a contribution of more than 5% to 
the total curvature (notice that this also comes with a 
change in the algebraic properties of the spacetime, since 
the principal null directions "peel apart" as each be- 
comes important with decreasing r). In other words, as 
the source is approached from r = cx), Z„ is the region 
in which starts to make a significant contribution. In 
the case of our BBH simulations, we find that ^'4 consti- 
tutes more than 95% of l^*] in the region r > 200M. We 
observe that ^"3, ^2 and ^'i begin to contribute > 5% to 
the curvature at r = 200M, 75M and 15M, respectively. 
Note that the precise values of r depend on the curve 
along which the falloff is measured and on the choice of 
cut-off percentage. 



FIG. 3. The fractional contribution of each (measured 
along Ai) to the total curvature, which we define as l^*] = 
X]t=o zones Zn are the regions in which '^'k, for 

k > n, give a contribution of more than 5% to Beyond r > 
200M, dominates and may be reliably used as a measure 
of gravitational radiation. Note that Zq does not appear as 
^/o is already below 5% at r = lOM, where we can first start 
tracking the peak in \'^4\. 



IV. DISCUSSION 

We have performed a 3-orbit BBH simulation and mea- 
sured the falloff of the Weyl scalars, obtaining results in 
agreement with the peeling theorem to within 4%. 

There are many approximations introduced in convert- 
ing the precise assumptions of the theorem into practical 
numerical calculations. For example, in this work we 
approximated null geodesies by tracking the location of 
peaks in 1^*4 1 along the x-axis. This neglects any angular 
component in the null geodesies and also assumes that 
the peaks propagate along null geodesies. Furthermore, 
we used a coordinate tetrad which we assume satisfies 
the assumptions of the peeling theorem. The fact that 
we found agreement with the predictions gives strong ev- 
idence in support of the approaches and approximations 
typically used in NR simulations. 

As in [T], accurate falloff rates were easily extracted 
for ^'2, ^3 and ^'4, where the finite differencing error 
was negligible and the deviation from the expected rate 
was almost certainly due to computing the falloff at a 
finite radius, where sub-leading terms are nonzero. 

However, due to the presence of a large amount of nu- 
merical noise, and ^'i proved much more difficult to 
analyse. This is not surprising given that the spacetime 
we are considering — a BBH inspiral — has a strong out- 
going radiation component {'^4 and "^3) and Coulomb- 
type potential (^'2) and only weak incoming radiation 
(^'i and 5*0) • Nevertheless, by minimising errors com- 
ing from inter-patch and mesh refinement boundaries, we 
were able to compute and in a sufficiently large 
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region to extract falloff rates within 4% of the expected 
values. 

As discussed in Sec. |III| and Fig. [2] the rapid falloff 
of vj/g means that it drops to a level where it is strongly 
affected by numerical noise for r ^ 60M. Furthermore, 
for r ^ lOOAf it reaches a point where it is orders of 
magnitude below the other '^^ and its contribution to 
the curvature is negligible. This indicates that any ap- 
proximations based on the vanishing of \l/o — for example 
"freezing-^'o" boundary conditions |54|, I55j — are robust 
at these large radii. 

We have also shown a representative example of the 
radial zones in which each of the Weyl scalars begins to 
contribute to the total curvature. In the case of our BBH 
simulations, we found that ^'4 constituted more than 95% 
of 1 5* I in the region r > 200A/. The identification of the 
zones in Fig. [3] is valid for the equal-mass, non-spinning 
BBH configuration we have studied here. While we ex- 
pect this to be representative of other configurations — 
possibly involving spins and unequal masses — it is likely 
that the exact locations of the peeling regions will be 



different. 

The use of ^4 to compute gravitational waveforms 
from NR BBH solutions is based on the assumption that 
the peeling theorem can be applied. Although it is rea- 
sonable to expect that this is true, this work provides 
reassuring confirmation that this is indeed the case. 
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